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Abstract: We report an optimal extraction methodology, for the reduction of multi-object fibre spec- 
troscopy data, operating in the regime of tightly packed (and hence significantly overlapping) fibre 
profiles. The routine minimises crosstalk between adjacent fibres and statistically weights the extrac- 
tion to reduce noise. As an example of the process we use simulations of the numerous modes of 
operation of the AAOmega fibre spectrograph and observational data from the SPIRAL Integral Field 
Unit at the Anglo-Australian Telescope. 
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1 Introduction 

Fibre-optic multi-object spectroscopy is a powerful tool 
with which to segment the focal plane of a wide-field 
spectroscopic teles c ope. The technique, pioneered in 
the 1970s (see iHilll (119881 1 for a historical review) , has 



been repeatedly put to use for major undertakings in 
astronomy. The high multiplex advantage that can 
be achieved, typically over wide fields-of-view, using 
fibre optic feeds easily outweighs the expense of im- 
plementation and inherent limitations of a fibre feed 
system for many classes of astronomical observation. 
Notable recent e xamples of in clude the 2dF Galaxy 
Red shift Survey (ICoUess et al\\2001i ) . 2dF Quasar sur- 
vey 
vey 



Bovle et alll2000l') and the Sloan Digital Sky Sr 
York et alhood ). 



The use of coherent fibre bundles to record two- 
dimensional spatially resolved information at the spec- 
trograph slit is an obvious extension of the technol- 
ogy, and indeed a number of the current generation 
of Integral Field Spectrograph (IFS) systems use fibre 
feeds to reformat the focal plane (e.g. AAT-SPIRAL, 
Gemini-GMOS, ESO-VIMOS & ARGUS). 

However, increased multiplex (or field of view for 
IFS instruments) comes at a price. With CCD (and 
IR array) pixels still at a premium in any astronom- 
ical instrument, there is a need to tightly pack the 
fibres into the spectrograph slit, minimising the dead 
space between adjacent fibres in order to maximise the 
number of fibres available for a given instrument for- 
mat. But tight fibre packing means poor sampling 
of individual fibre Point Spread Functions (PSF) and 
can result in significant overlap of the profiles of ad- 
jacent fibres leading to strong cross-contamination, or 
cross-talk, between fibres. Recovering the spectral in- 
formation for spectra observed in tightly packed sys- 
tems is therefore complex. Simple summation of the 
pixel values surrounding the peak of a fibres' trace is 
doomed to failure due to contamination from adjacent 
fibres (particularly with high contrast observations as 
are common in IFS applications). Even if individual 
fibres are reasonably well resolved, the extraction will 
often suffer increased noise from contamination by the 



wings of adjacent fibre profiles, or poor weighting of 
pixel values introducing increased CCD read-noise in 
the limit of low background observations. 



What is required is an extraction procedure which 
takes account of the interaction between adjacent fi- 
bres and provides a statistically optimal estimate of 
the true intensity for each fibre spectrum. We have 
recently implemented such an extraction mechanism 
for the AAOmega spectrograph system at the Anglo- 
Australia Telescope. The code has been developed 
in the IDL programing language and implementation 
within the 2dfdr data reduction environment (con- 
structed primarily in FORTRAN). In this work we quan- 
tify the need for such an extraction algorithm and 
demonstrate that while there is limited gain in us- 
ing the routine for the well resolved fibre profiles of 
the default AAOmega instrument Multi-Object Spec- 
troscopy (MOS) mode, it is essential for high purity 
data taken with the tight fibre packing utilised for the 
recently implemented AAOmega mini-shuffling nod- 
and-shujfle observing mode and for the SPIRAL IFS 
feed to the AAOmega spectrograph. We describe the 
methodology used and demonstrate its application to 
observations taken with the AAOmega-SPIRAL sys- 
tem. 



The paper is organised as follows. In § [2] & [3] 
we introduce the three alternative spectra extraction 
methodologies currently implemented with in the 2df dr 
data reduction software and four modes of instrument 
operation for the AAOmega spectrograph system at 
the A AT. In §|4]we demonstrate the need for an extrac- 
tion algorithm which accounts for fibre-to-fibre cross- 
talk. The two non-trivial extraction algorithm from 
§ [2] are then described in detail in § [S] The proce- 
dure for applying the the optimal extraction algorithm 
is given in §[n]&[Zl Scattered light effects are addressed 
in § [8] and example data from the AAOmega-SPIRAL 
system is presented § ^ 
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Spectral extraction method- 
ologies 



Three spectral extraction methodologies are currently 
implemented within the 2df dr data reduction environ- 
ment used with data from the AAOmega spectrograph. 
In each case the starting point is a map of the fibre 
profile centroids on the CCD. Since this fibre map is 
a series of near parallel tracks it has historically been 
labeled the tramline map. The generation of the tram- 
line map is conducted separately from the spectral ex- 
traction and is not discussed in detail here. Briefly, 
the fibre profile centroids are determined at intervals 
across the CCD (every ~50 pixels along the dispersion 
axis) and a low order polynomial model is fitted to the 
centroids for each fibre, guided by an optical model for 
the expected spectrograph camera distortions. 

The three spectral extraction methods currently 
available within 2dfdr are: 

Tramline summation - This most basic of extractions 
is obtained by simple summation of all pixel values 
associated with a given fibre. The spectra for the fi- 
bres are dispersed broadly along rows of CCD pixels 
and each column of the CCD is treated independently. 
The summation range for each fibre typically runs over 
the pixels bounded by the mid-points between the two 
adjacent fibre profiles. While quick to compute, the 
simple tramline extraction propagates the maximum 
CCD readout-noise into the final extraction since it 
gives equal weight to all pixels in the summation re- 
gardless of the fiux level of the fibre profile in a given 
pixel. Tramline summation also sufi^ers an aperture 
loss effect if the inter-fibre gap over which pixels are 
summed is not significantly larger than the fibre pro- 
file width. This aperture correction can of course be 
accounted for if the true fibre profile is known. 

Gaussian weighted summation - An algorithm for per- 
forming a weighted Gaussian summation over a sin- 
gle isolated fibre profile (via a Least Squares Fit) is 
discussed in section 15.11 This weighted summation 
minimises the contribution for CCD readout noise and 
does not sufi'er an aperture effect. This mode is the 
default for AAOmega-MOS spectroscopy and we show 
in ^that for well separated spectra it delivers accept- 
able results. 

Multt-fibre deconvolution - The ideal solution for the 
minimisation of cross-talk between fibres is an extrac- 
tion algorithm which performs a multi-fibre deconvolu- 
tion of the data given an underlying model assumption 
for the fibre profiles on the CCD. 

The extent to which fibre-to-fibre cross talk is present 
under each extraction models described above is tested 
in ^ The algorithms adopted for the Gaussian and 
multi-fibre deconvolution solutions are then discussed 
in HS.ll and i)5.2l All three approaches have been imple- 
mented for the AAOmega spectrograph at the Anglo- 
Australia Telescope. The multiple operating modes of 
AAOmega are outlined in the next section. 



3 The fibre spectrograph in- 
strument modes 

Table[T]outlines four instrument modes fo r the AAOmega 
facility at the Anglo- Au stralia Telescope (jSaunders etali 
l2004l : (sharp aZfeoO^ I . Data from the AAOmega sys- 
tem is primarily processed using the 2df dr data reduc- 
tion software and the three extraction algorithms pre- 
sented above have been developed from the code base 
of the previous 2dF spectrographs. 



AAOmega MOS 

The first instrument mode is the default multi-object 
spectroscopy (MOS) mode of operation. In this mode 
392 science fibres are deployed on astronomical targets 
across the vrdeg^ field-of-view of the 2dF prime focus 
corrector. The fibres feed the dual beam AAOmega 
spectrograph, each arm of which is equipped with a 
2kx4k E2V CCD. This resuhs in a fibre-to-fibre sep- 
aration, the fibre pitch, of ~10 pixels. The cameras 
deliver largely Gaussian profiles of FWHM~3.4 CCD 
pixels. 



Mini-Shuffling 

The second mode under consideration is the newly 
implemented mini- shuffling mode. The high quality 
AAOmega PSF makes it possible to undertake fibre 
nod- and- shuffle observations for high quality s ky sub- 
traction (|Glazebrook fc Bland-Hawthornll200ll ) by in- 
terleaving multiple on-sky exposures on a single CCD 
frame using the inter-fibre gaps to provide the required 
CCD storage areas during charge shuffling. This mode 
effectively doubles the number of AAOmega fibres on 
a CCD by halving the fibre pitch. This introduces sig- 
nificant fibre-to-fibre crosstalk as will be demonstrated 
in 21 A subsequent paper will discuss the effectiveness 
of mini-shuffling in detail. 



The SPIRAL Integral Field Unit 

The optimal extraction algorithm presented in i)5.2l 
was developed for the AAOmega-SPIRAL Integral Field 
Spectrograph system used with AAOmega at the AAT. 
It is based on an impleme ntation developed for the 
CIRPASS IR spectrograph (|Parrv efa/.l 120041 '). Both 
CIRPASS and SPIRAL use a tight fibre packing, and 
with good reason. CIRPASS initially utilised a Rock- 
well Hawaii-IK IR array and hence required a tight 
fibre p acking to maintain the wide field o f view in IFS 
mode l|Krainovic. Sharp fc Thattd 120061 ) . The SPL 



RAL system is designed to allow nod- and- shuffle (jGlazebrook fc Bla: 
bOOlh observations with the full 22.4x 11.2 arcsec FoV 
of the SPIRAL system, while using the same detec- 
tor real-estate as the AAOmega multi-object fibre feed 
from the 2dF fibre positioner (Sharp et al. 2006). Hence 
for the SPIRAL system CCD pixels are at a premium 
as each fibre essentially requires twice as many pixels 
on the CCD as dictated by its intrinsic foot print. 
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HERMES 

The final mode of operation considered in Table [T] is 
that proposed for the HERMES high resolution (R~30,000) 
multi-object stellar spectrograph under development 
for use with the 2dF positioner at the AAT. HERMES 
will employ 4kx4k CCDs with 15/im pixels. However 
at /1. 5, the HERMES camera design is slower than 
the /1. 3 cameras of AAOmega, resulting in a broader 
PSF (5 pixel) than that of the currently implemented 
AAOmega modes. This results in an increased profile 
overlap compared to the basic AAOmega mode. 



Table 1: The pertinent parameters for the four 
model instrument modes of the AAOmega system 
are given. 



Instrument 


Fibres 


Pitch 


FWHM 


Mode 




(pixels) 


(pixels) 


AAOmega 


392 


10.0 


3.4 


Mini-Shuffle 


784 


5.0 


3.4 


SPIRAL 


512 


4.0 


2.4 


HERMES 


392 


10.2 


5 



4 Fibre-to-Fibre cross contam- 
ination 

Before presenting our new extraction algorithm which 
reduces fibre-to-fibre spectral cross-contamination (5 15. 2|) 
we will motivate its' development by consider the level 
of this cross-talk in models of the four observing sys- 
tems of Table [1] The pertinent parameters for the 
models of each observing mode are: 

Fibre pitch - The spacing between fibres on the CCD. 

Fibre profile - e.g. the projected FWHM of a fibre pro- 
file on the CCD. 

Relative intensity - The intensity ratio between spectra 
to be extracted. 

Fig.[T]presents a visual representation of the instru- 
ment modes of Table [1] A pair of isolated fibres are 
considered. This represents the limiting case of cross- 
talk between two adjacent fibres only. Two model fibre 
profiles are shown with an input A mag=3, a fiux ra- 
tio of ratio~15.8. The fibre centre, and ±3(t range 
are marked for each fibre (solid bars) as are the inter- 
fibre ranged associated with each fibre (dashed bars). 
This later range is used for the Tramline summation 
extraction, while ±3(t is used for the Gaussian sum- 
mation unless it is larger than the inter-fibre range, 
in which case it is truncated to this range. A visual 
inspection of Fig. [1] indicates that for all but the de- 
fault AAOmega mode there may be large fibre-to-fibre 
cross-contamination, due to the significant overlap be- 
tween fibre profiles, if a simple pixel intensity summa- 
tion extraction is used. For the later instrument modes 
an extraction process is needed that accounts for the 
fiux of both fibres simultaneously and attempts to cor- 
rectly distribute the intensity information between the 
two fibre profiles in the overlap region. A quantitative 
analysis follows. 

Examination of observational data from the AAOmega 
system indicates that a Gaussian fibre profile is an 
acceptable first approximation. We will use the Full 
Width at Half Maximum, FWH\£|, as the measure 
of profile width in what follows. We express relative 
spectral intensities in terms of the magnitude differ- 
ence of the integrated profile fiuxes, A mag, for two 

^We define the inter-fibre range for each fibre as ithalf 
the fibre pitch, centred on the fibre. 

^The FWHM is related to the Gaussian width, a, as 
FWHM=2v'2 In 2 X ct ~ 2.355a 



astronomical objects illuminating the fibres. Conven- 
tional wisdom for the AAOmega system has been to 
limit Amag<3. 

Fig.Eldemonstrates that for the default AAOmega- 
MOS mode of operation (the two top figures in each 
column) cross-contamination of the fainter spectrum 
is limited to <1.6% of its integrated profile flux for 
Amag<3. Improved accuracy requires a reduced A mag. 
However the three remaining instrument modes de- 
tailed in § |3] all suffer significant cross-contamination. 
In these cases not only is cross-contamination of the 
weaker spectrum increased significantly, but double 
counting of flux leads to >10% errors in fiux normalisa- 
tion when A mag<0.3. Clearly an improved extraction 
methodology is needed which can account for the over- 
lapping fibre profiles in these cases. No figure is given 
for the Multi-profile deconvolutton extraction since, un- 
der the assumption of this idealised initial analysis, it 
returns perfect results for all four instrument modes. 
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Figure 1: Model fibre profiles are shown for a pair 
of fibres with each of the instrument modes given 
in Table [TJ The logarithmic scaling heightens the 
visibility of the region of overlapping profiles. The 
profile intensity ratio is ^^15. 8 giving Amag=3. 
The fibre centres, and ±3a range are marked for 
each fibre (solid bars) as are the free inter-fibre 
ranges associated with each fibre (dashed bars). 
For all but the first instrument profile, the 3(t 
range extends beyond the inter-fibre range. 
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Figure 2: Using two model Gaussian fibre profiles and the instrument modes of Table [T] the ratio of 
extracted flux to input model profile flux is presented as a function of the relative intensity difference 
(expressed as a magnitude difference, A mag) between the two input profiles. The left column shows 
results using the Tramline summation extraction, the right the Gaussian least squares fit. In all cases the 
second fibre is A mag fainter than the brighter first fibre. An accurate extraction would follow the 1:1 
locus. For all instrument modes the extracted flux for fibre 2, the fainter profile, is increasingly in error 
as A mag increases due to cross-contamination from the brighter source fibre. No traces are shown for 

the multi-profile deconvolution extraction as the process is found to be 100% accurate (the correct flux is 
recovered for both fibres) in the noise free model limit for all four instrument modes. 
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4.1 Effects of errors in centroids and 
profile widths 

Two parameters of the assumed fibre profiles domi- 
nate the extraction error, errors in the assumed fibre 
centroid and errors in the profile widths. Fig.[3]demon- 
strates the effects of introducing random modifications 
to the fibre centroids and widths (in isolation). Fifty 
thousand realisations are computed with a random er- 
ror introduced to each fibre profile drawn from a nor- 
mal distribution with the indicated width. The ratio 
of the mean extracted fiux to the input model value 
(and associated il-cr scatter) is given for each error 
distribution. The fiux ratio between the fibres was 
set to unity (Amag=0). The simulations show that 
while errors in the fibre centroid value may be toler- 
ated at the 0.2 pixel level (6% of the FWHM here), 
an accurate knowledge (<0.1 pixel, ~3% of FWHM) 
of the fibre profile is required to control errors under 
all extraction methodologies. 

4.2 Effect of Signal-to-Noise ratio 

A final parameter to consider when testing the accu- 
racy of the extraction of fibre profiles is the signal- 
to-noise ratio (SN) of the data. To investigate this, 
he two fibre profile model is retained, but a constant 
readout noise per pixel and a shot-noise component 
is added to each realisation of the model data based 
on the intrinsic strengths of the input model profiles. 
Profiles are generated with a range of normalisations 
while the input magnitude ratio (A mag) is held fixed 
for each realisation. Fig. |4] shows the results obtained 
for the three extraction methodologies, for a range of 
SN ratios, assuming the default AAOmega-MOS mode 
as described in Table [1] 

The figure shows the intrinsic signal-to-noise ra- 
tio used when generating the model profiles and the 
signal-to-noise ratio of the extracted profile flux esti- 
mate. This recovered signal-to-noise is deflned as the 
ratio of the extracted count rate to the noise level 
as estimated from the model data by the extraction 
process. When reviewing the results of the simula- 
tions, one should note that the SN ratio recovered for 
each spectrum is a statistical measure derived from the 
data, and does not directly address the issue of the 
double counting of fiux due to fibre-to-fibre crosstalk. 

Fig. 13] indicates that for the default AAOmega- 
MOS mode the choice of extraction methodology has 
only limited effect on the final SN ratio. In the case 
of Amag=0 there is some indication that the multi- 
component decomposition (optimal extraction) reduces 
the scatter in the ratio between recovered and intrinsic 
profile flux. This increase becomes more pronounced 
for lower SN ratio data. The situation is similar for 
Amag=3 with the error in the flux value derived for 
the fainter fibre suffering a greater inflation for the 
tramline and Gaussian extraction in comparison to the 
optimal extraction. 

The process is repeated for each of the instrument 
configurations from Table [T] Fig. [S] shows the results 
for Amag=0 and Fig. [6] for Amag=3. It is immedi- 
ately apparent that in the case of tightly packed fibres 
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Figure 3: The ratio of the mean measured flux 
to the input model intensity (and associated RMS 
scatter) is shown for the three extraction meth- 
ods of § [5] after 5,000 independent reahsations of 
the error model. For each realisation the model 
fibre profile centroid and FWHM parameters are 
drawn from normal distributions with the indi- 
cated widths. The summation aperture correction 
factor has not been applied to the Tramline ex- 
traction, leading to its' offset below the other two 
extraction methods. 
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flbre-to-fibre cross-talk plays a significant role in mod- 
ifying the extracted flux values for the tramline and 
Gaussian extraction methods. 
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Figure 4: Variations in the extracted flux are considered as a function of the signal-to-noise ratio intrinsic 
to the model date realisation. Signal-to-noise is quoted for the integrated profile flux. Two sets of 
simulations are shown, both for the AAOmega-MOS instrument mode. Fifty thousand independent trials 
of the noise model are run with the fibre profile fiux ratio set to Amag=0 & 3 (top and bottom rows 
respectively). Fibre 2 is the fainter fibre in the Amag=3 case. The left column presents the ratio of the 
extracted flux to that intrinsic to the model. The right column compares the intrinsic SN to that recovered 
for the extracted profiles. In the case of A mag=0 the results for the two fibres are indistinguishable. The 
optimal extraction is seen to return a smaller scatter in the recovered flux ratio, the scatter remaining 
largely constant as SN decreases. For the Amag=3 case both the Tramline and Gaussian extraction 
return an increased scatter in the ratio for the fainter fibre (fibre 2) of the pair. The RMS error bars for 
fibre two in the lower right plot are similar for all three extraction algorithms and so for clarity only those 
associated with the OptEx algorithm are given. 
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Figure 5: Fig. 0] is repeated for each instrument configuration presented in Table [T] and with Amag=0. 
For the compact fibre systems a significant systematic flux error is seen (for both fibres) in the extracted 
flux due to double counting in the overlap regions by the Tramline and Gaussian extraction methods. 
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Figure 6: Fig. [5] is repeated with Amag=3. Even with a large flux ratio between the fibres the Optimal 
extraction method recovers the correct flux for both the bright (fibre 1) and faint (fibre 2) profiles. The 
apparent increase in the SN ratio of the faint fibre with the tramline and Gaussian extraction is due to 
the erroneous double counting of flux from the adjacent flbres, which would represent a degradation of 
the SN ratio actually achieved. 



www.publish.csiro.au/journals/pasa 



11 



5 Profile extraction algorithms 

In each of the three spectral extraction methodologies 
introduced in 321 it is assumed that each fibre projects 
a spectrum onto the detector such that the light is dis- 
persed along rows of pixels. At each spectral pixel the 
fibre profile is assumed to be well modelled by an an- 
alytic function aligned precisely along columns of the 
detector, orthogonal to the dispersion axis. These as- 
sumptions are approximately true at the pixel-to-pixel 
level, the macroscopic spectral curvature on the CCD 
being on a much larger scale. Such an approximation 
is made in almost all approaches to the extraction of 
spectra from 2D CCD data. 

In the discussion of the two non-trivial spectral ex- 
traction algorithms presented below we use the terms 
defined in Table H 

5.1 Gaussian summation extraction 
via least squares 

An extraction method which minimises the read-noise 
contribution from pixels in the wings of the fibre pro- 
files, while maintaining a full summation of all of the 
observed fiux, is in principle rather trivial for the case 
of a single simple spectral profile (although careful con- 
sideration of the system is required for a truly optimal 
result, see iHornd (|l986t ) for a comprehensive discus- 
sion) . 

We consider each column of the dispersed data in 
isolation, and construct a model for the intensity at a 
given spatial pixel assuming an appropriate fibre pro- 
file (/f>i, and integrated intensity i]. 



Mi ^T) X 



(1) 



In the case of non-overlapping spectra, and given a 
recorded CCD pixel intensity Di with associated error 
estimate (jJ3, an estimate of the integrated profile 77 is 
given by the minimisation of 



R 



E 



This yields 

dR _ sr^ 
duk ^ 

i 

and on minimisation, 



■tih 



E 



2nk(f>i 



(2) 



(3) 



(4) 



with the summation running over a sufficiently large 
range of pixels, i, to cover the full fibre profile, without 
overlapping with adjacent profiles. A range of ±3(t 



''The error estimate for each data value is, as is typical 
for such data, constructed from the data under assumptions 
of the detector read-noise and gain characteristics and the 
shot-noise from each data value itself. While this clearly 
leads to an imperfect estimate of the error in the observa- 
tion, we find this simple approximation to be sufficient for 
our requirements. 



(or half the inter-fibre separation for tightly packed 
fibres) is typically adopted. The resulting value of the 
integrated profile fiux rj is calculated for each fibre in 
isolation and for each element along the dispersion axis 
in tern. 

5.2 Multi-profile deconvolution ex- 
traction 

The core of the fibre profile achieved by the AAOmega 
system is observed to be well matched to a Gaussian 
profile; and for speed and simplicity, while still main- 
taining sufficient accuracy, the current 2df dr data re- 
duction environment used with AAOmega data imple- 
ments a simple Gaussian profile. 

Furthermore, we assume that in the first instance 
the precise position of each fibre profile as a function 
of wavelength (i.e. detector column) is known from the 
previous trace of an appropriately illuminated fiat field 
frame as discussed in ^ 

We then define a normalised spatial profile for fibre 
k, as a function of the CCD pixel i along the spatial 
axis (the fibre slit axis, orthogonal to the dispersion 
axis), such that 

i 

Hence the contribution of fibre k to the observed 
intensity value at CCD pixel i is given by the product 
of the integrated profile intensity for the fibre, r/k, and 
the fibre profile at pixel i. For any given pixel of the 
detector a model value for the recorded count rate at 
pixel i is then simply the sum over all contributing 
fibre profiles, 



(6) 



Given Di and at, the recorded count rate and an 
estimate of the statistical error at each pixel (along a 
single column of the detector), we wish to evaluate rjk 
for each spectrum. This is achieved via the minimisa- 
tion of the residual, 



^ 1 (A - M,y 

^-2Z. ^ 



(7) 



and on setting — and substituting for Mi we 
find. 



E\ ^ 4>ji4>ki 



zA 



(8) 



Letting Ckj=Y,, -7^^ and hj=Y.i 
we find 



D^4>Ji 

„2 



E 



(9) 



Equ.[5]can be solved using any of the multitude of 
methods for solving coupled linear equations. Numeri- 
cal stability considerations have guided our implemen- 
tation towards a Singular Value Decomposition (SVD) 
approach. SVD allows terms of low significance to be 
removed from the matrices during fitting, providing 
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Table 2: Definition of terms used in the discussion of tfie extraction algorithms. The desired products of 
the fitting process are the integrated fibre intensities rjk for each fibre k, and the associated variance in 
T]k expressed as vark- The process includes all pixels, i, along the spatial (slit) axis of the CCD which is 
assumed orthogonal to the dispersion axis. It is then repeated for each pixel of the dispersion axis. 



Parameter 


Definition 


fib 

i 

k 
j 


The number of fibre profiles in the system 

Running index over pixels in the spatial (CCD y-axis) 

dimension orthogonal to the dispersion axis 

Running index over successive fibre profile l-Nfu, 

Dummy index for fibre k when computing profile overlaps 


Di 
M, 

Var{D,) 


CCD data value at pixel i 
Model for data at pixel i 

For small fibre pitch this is a sum over overlapping profiles 
Error estimate at CCD pixel i 

Variance estimate at CCD pixel i with data value Di 


Vk 

vark 
't'i 

4>ki 


Integrated intensity for fibre k at the current wavelength element 
Computed estimate of the variance associated with fibre k with fiux i]k 
Normalised fibre profile for a single fibre at spatial axis pixel i 
Normalised fibre profile of fibre k at pixel i 


^kj,c'f.j 


Internal variables for the cross terms between the profiles of fibre k 
and fibre j for all fibres k,j=l-N fn. 

Internal variables for the sum of products of the fibre profile j 
and the data/ variance values at pixel i 



control over numerical stability for high contrast data 
sets. One notes that, even for very tightly packed fi- 
bre data, the overlap between fibres separated by more 
than a few FWHM will be close to zero in most in- 
stances (with the exception of very high contrast ratio 
data, or for a model PSF with wide scattering wings). 
This makes c^j at least band diagonal, and tri-diagonal 
in many cases, allowing trivial solutions to Equ. [9] in 
these cases. 

5.2.1 Propagation of variance information 

Propagation of error information is critical for all ob- 
servational science, but is woefully missing in many 
spectral extraction routines. It can be achieved as fol- 
lows. For each pixel i in the input image we calculate 
the total of all fibre profiles, (j)ki, that contribute to 
that pixel, 

^^=^r)k<t>k^ (10) 
k 

Then for each spectrum k, we calculate the fractional 
contribution of this spectrum's profile to pixel i, 

F.. = ^ (11) 

An estimate of the total variance in the extracted in- 
tensity r]k is then given by 

varfc=^FLxa? (12) 



5.2.2 Alternate variance solution 

The variance estimate associated with any given de- 
tector pixel, Var{Di) is determined directly from the 
observed pixel value and the predetermined read-noise 
and gain parameters for the CCD system. Given Equ. [6] 
the value is given by the sum of the contributions from 
each fibre at that pixel plus the additional detector 
read-noise. Under the assumption of shot-noise in the 
photon arrival rate at a pixel, the contribution of a 
particular fibre to the variance estimate for that pixel 
is given by the produce of the total count rate for the 
fibre and the normalised profile intensity at the pixel. 

Vi = 4>kir)k = (t)kivark (13) 

When multiple fibre profiles contribute to a pixel, a 
model for the observed variance in that pixel is then 

= (j}k^var}j + (14) 

We then wish to minimise the residual between this 
model and the variance in the observed data, Var(_Di). 

i 

(16) 

By analogy with Section 15.21 and solving or for 
(5R/5iiarfc=0, we therefore arrive at an estimate for the 
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variance, vart, associated with each integrated fibre 
profile intensity rjk by solving 

^war-fcc'fcj = b'j, (17) 

k 

with ci~J2i i't'kift'ji) and b'j=J2i {Var{Di) ~ n^^) 
This can be solved directly alongside Equ.|9] with min- 
imal additional computational overhead. 



5.2.3 Practical limitations and the itera- 
tive solution 

As demonstrated in Fig. O variation in the fibre pro- 
file FWHM with wavelength is of particular concerr0. 
Poor correction for background scattered light will also 
prevent a satisfactory solution being achieved. Results 
can be significantly improved via an iterative solution 
to the problem of determining these free parameters. 

We adopted an iterative solution with a number of 
free external fitting parameters. The model is derived 
from a high signal-to-noise flat field frame and then 
the parameters are locked at the resulting values when 
fitting for science frame intensities. The free parame- 
ters are 

Background - We include a low order polynomial fit to 
the broadband scattered light distribution across the 
CCD (pedestal plus gradient). In this model, scattered 
light is effectively treated as a pedestal correction to 
the observed data. For the AAOmega system, the scat- 
tering is well modelled as a DC pedestal offset across 
the CCD. This background level is a free parameter 
in the iterative solution of Equ. [5] (but see § [7] for a 
further refinement). 

Fibre Profile - The normalised fibre profile under the 
assumption of a Gaussian PSF has only a single free 
parameter, a, the Gaussian profile width. A more ac- 
curate accounting for the SPIRAL fibre PSF is consid- 
ered in § 1121 

2D parameter distribution - A full 2D polynomial fit to 
the variation of these free parameters as a function of 
CCD position should also be performed. This can be 
achieved by fitting a subset of data columns at regular 
positions across the dispersion axis. A low order model 
would then be fitted for the full data set guided by 
instrument design considerations. This has not yet 
been implemented within the 2dfdr software. 

With these external parameters in place we iter- 
atively solve Equ. for all fibres simultaneously, for 
each element of the dispersion axis in turn. T he IDL 
prototyping code used MPFIllf] ( Markwardtll2009l '). adopt- 
ing the statistic as the figure of merit applied to the 
difference between the data and model. 



■'The AAOmega red camera suffers from a degraded 
(■-^O. 2 pixels) focus at long wavelengths which is attributed 
to the fast {f/1.3) camera and increased photon penetration 
into the CCD before detection. 

^http: / /cow. physics, wise. edu/^craigm/idl/fitting. html, 
MPFIT, an excellent IDL implementation of the Levenberg- 
Marquardt fitting procedure (iMarkwardtii2009l) 



5.2.4 Run time 

The multi-fibre deconvolution algorithm described above 
is computationally intensive and to be of value the 
time taken to process observational data must be con- 
sidered. The AAOmega system (Table [l| requires the 
solution for ~400 simultaneous fibre profiles (~800 for 
the mini-shufSe mode) across 4000 spatial pixels for 
each of the 2000 spectral columns of a 2kx4k E2V 
CCD. The algorithm has been successfully implemented 
in the 2df dr software environment (written primarily 
in (Fortran95) and the extraction is completed for a 
single frame in ^l-2minutes on a modern desktop PC. 

In its current format, fifty percent of the calcu- 
lation time for each instance of the fitting process is 
concerned with the calculation of fibre profiles and 
the population of the Cjk & bj matrices from Equ. [S] 
The current implementation does not distribute calcu- 
lations across multiple PCUs and so this minor code 
extension should allow for a near linear speed increase, 
for these calculations, proportional to the a number of 
processors available on a modern multi-core computer. 

6 The procedure 

For data taken with the SPIRAL- AAOmega system we 
undertake the following procedure 

• A fibre-flat-field frame is taken in which all the 
fibres are illuminated with a uniform continuum 
source in order to trace the centroid of the fibre 
profiles across the CCD. Fibre-fiat-field frames 
can also be used to correct for variations be- 
tween fibres in the relative response function 
with wavelength. They do not provide an accu- 
rate correction for pixel-to-pixel CCD response 
variations due to the high degree of spatial struc- 
ture across the fibre profiles. If such a correction 
is required, e.g. to correct for interference fring- 
ing in the CCD, a more uniform long-slit-flat 
field frame would be required. 

• All CCD frames are processed for overscan/bias 
correction and population of the variance array 
information. 

• An iterative fitting algorithm is applied to ex- 
tract the integrated fibre profile intensities for 
the flat field. Experience with SPIRAL data shows 
that re-adjustment of the fibre centroids mea- 
sured from the fibre-flat-field is not required for 
AAOmega. 

• The iterative fitting algorithm is applied to the 
science data, holding the centroid and PSF val- 
ues fixed (determined from the fibre-fiat-field above) 
and fitting only for the pedestal offset as a free 
external parameter when iteratively solving for 
the integrated fibre intensities via Equ[51 Using 
the alternative background subtraction method 
outlined in section SJZI below, the extraction of 
science data can be made non-iterative once the 
extraction parameters have been determined from 
the fibre-flat-field data. 



14 



Publications of tlie Astronomical Society of Australia 



7 Non-iterative background so- 
lution 

If a simple pedestal/polynomial model is adopted for 
the background component then it is possible to fit the 
background without resorting to an iterative solution. 
For example, a simple pedestal can be thought of as an 
additional fibre profile with a uniform profile intensity 
at each pixel. 

Higher order polynomial terms of the form could 
also be introduced, up to some number of terms Nbg- 
One would then solve Equ.[9]for rjk with 1 < fc < Nfib 
for each fibre with normalised fibre profile (f)ki and also 
for Nfu, + 1 < fc < Nfu, + Nbg governing the back- 
ground model intensity terms. 

This more elegant formalism allows the solution 
of Equ. [9] without iteration in the case that all other 
external variables are predefined. 



8 Scattered light 

Two basic assumptions of the spectral extraction pro- 
cesses discussed in this work are: 

1) The spatial (slit) axis on the CCD can be considered 
orthogonal to the dispersion (spectral) axis at each 
pixel along the dispersion axis 

2) Each spectral pixel is treated as if it where indepen- 
dent of the adjacent spectral pixels, i.e. no accounting 
is take of the 2D nature of the instrument PSF. 

Essentially we are replacing the real- world spectro- 
graph Point Spread Function (PSF) with a pixel-by- 
pixel Line Spread Function perpendicular to the dis- 
persion axis. 

Such assumptions are made primarily for compu- 
tational expedience. For a compact and well-sampled 
PSF they are valid. This is fortunate since a full 2D 
modeling, which amounts to simultaneous deconvolu- 
tion involving every pixel on the CCD (and a number 
of virtual pixels extending beyond the physical device), 
will likely remain beyond the computing power avail- 
able for routine multi-fibre observations for some time 
due to the inherent scale of such a problem. As an 
example of this scale consider that the 512 fibres of 
the AAOmega-SPIRAL system each have ~2000 spec- 
tral elements (~830 independent resolution elements) 
covering ~8 million CCD pixels. 

The model breaks down if the fibre PSF possesses a 
measurable scattering wing which extends significantly 
beyond the core of the fibre profiles. If such a broad 
scattering profile is present then the assumption that 
data can be accurately modeled as a ID line spread 
function rather than the 2D PSF breaks down since 
scattered light present at any given position along the 
dispersion axis will be largely dominated by the inten- 
sity of the local average spectral intensity of the sur- 
rounding spectrum and unrelated to the specific spec- 
tral intensity of the current spectral pixel. 

In high signal-to-noise data (including the fibre- 
flat-field frames for science data that will ultimately 



be in the low signal-to-noise regime), even a scatter- 
ing fibre profile with low overall percentage light level 
in the scattering component can have significant side 
effects on the ultimate data quality if the profile has 
a broad PSF. The scattered light from an increasing 
number of more remote fibre profiles will become more 
and more important. 

8.1 Scattering without fibre struc- 
ture 

Some scattered light is inevitable in all optical sys- 
tems. Low level background structure which is largely 
devoid of the high spatial frequency signature of the fi- 
bres themselves can be removed relatively easily. One 
may fit directly to the background level observed in 
regions free from contamination by the fibre profiles. 
Alternatively a low order model can be included as an 
additional fibre profile as demonstrated in [JT] This 
model profile would ideally be informed by an inde- 
pendent modeling of the background scattered light 
component. 

8.2 Multi-component profile model 

For scattering more directly associated with the ex- 
tended wings of the individual fibre PSFs, a more com- 
plex consideration is required. As previously men- 
tioned, a full 2D deconvolution across the entire CCD 
array is likely beyond the capabilities of modern sys- 
tems, at least at a rate comparable with the data rate 
from modern multi-fibre spectrographs. 

A solution is to maintain the ID orthogonal extrac- 
tion assumption, but using fibre profiles with extended 
wings. This approach is flawed however. For an ex- 
treme example of why consider the case of a strong stel- 
lar absorption line within a strong continuum source. 
The scattered light at line centre would be dominated 
by the strong local continuum average, due to the un- 
derlying 2D nature of the scattered light PSF. The 
absorption line core proflle would however be rather 
weak in comparison. This local variation in the spec- 
tral shape essentially induces a local variation in the 
assumed line spread function. One cannot merely fit 
for a core profile with a scattering wing modeled as a 
fixed percentage of the core profile flux. 

A solution is to model each fibre profile as a com- 
posite of two normalised profiles, one for the core pro- 
file, (fii, and a second broader scattering profile 0^. 
Each component will have an independent intensity, 
T)'^ and rj" . The model for the CCD pixel intensity at 
pixel i is then the sum, over all fibre profiles k which 
contribute to the pixel, of both the core and scattering 
profiles. 

Mi= ini^k^ + nUU, (18) 

Equ. 1^ is now solved for 2xNfHy fibre profile intensities 

There are inherent degeneracies in this composite 
profile description for each fibre. At low signal levels 
one can encounter significant problems with numerical 
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stability when solving for values of rik in real data. We 
have found the best results are achieved via a Single 
Value Decomposition (SVD) approach. The opportu- 
nity to suppress poorly behaved elements of the in- 
termediate solution matrices within the SVD solution 
(akin to discarding from the analysis the coupled linear 
equations which contain lim ited information, (as dis- 
cussed by iPress etal\ (dUi)) recovers numerical sta- 
bility and more than compensates for the increased 
computational burden presented by SVD over simpler 
solvers. The SPIRAL IFU data extraction shown in 
Fig. [7] were achieved using a composite profile extrac- 
tion and an SVD solver. 

8.3 Scattering from beyond the free 
spectral range 

An interesting side-efi^ect is encountered as one ap- 
proaches the edges of the CCD. For a spectrograph 
with a fiat response function in wavelength, scattering 
into an individual fibre profile from shorter and longer 
wavelengths will be much the same as that scattered 
out of the element in question (provided the spectral 
shape does not possess strong spectral breaks). How- 
ever, as one approaches the limit of the free spectral 
range of the spectrograph this symmetry no longer 
holds. The exigencies of spectrograph design will typ- 
ically place sharp spectral features such as those from 
order sorting filters, dichroic beam splitters or the grat- 
ing blaze just beyond the operational range of the spec- 
trograph. This means that the scattered light arriving 
in pixels towards the edge of the detector array may 
manifest a significant reduction in intensity since lit- 
tle light is scattered back onto the detector from just 
beyond the observed wavelength range. This has been 
found to be particularly true for the AAOmega system 
since the Volume Phase Holographic (VPH) gratings 
it employs have strong gradients in the blaze profile at 
the extreme ends of the intended observational ranges. 

9 Discussion 

A striking comparison of fibre-to-fibre crosstalk using 
the tramline, Gaussian and optimal extraction tech- 
niques is given in Fig. [T] 

The left panel of Fig. [7] shows the reconstructed 
image obtained using the simple tramline extraction 
model. The image is created by collapsing a 3D SPI- 
RAL data cube over the spectral dimension, to create 
a low resolution (0.7arcsec pixels) image of a standard 
star observation. The vertical structure is not a tele- 
scope diffraction spike, but rather the manifestation 
of fibre-to-fibre crosstalk. The SPIRAL fibres are ar- 
ranged in 16 banks (IFU short axis) of 32 fibres (IFU 
long axis). Hence all fibres in the central column of 
the IFU are subject to scattered light in the wings of 
the profiles of the fibres close to the bright star under 
observation. With a simple tramline extraction, which 
fails to account for this crosstalk, the central columns 
of the IFU exhibit an enhanced light level in all spectra 
in the central columns. For more exotic astronomical 
objects, any and all columns of the IFU with high con- 



trast data (either in the form of intensity or spectral 
type variations) will suffer a similar fate. 

For the central pane of Fig. [7] the extraction was 
performed using the Gaussian weighting. While this 
methods optimises the signal-to-noise ratio for each 
spectrum, it fails to account for the cross-contamination 
between the spectra. 

The right most image shows the same data after 
applying the optimal extraction technique. A round 
stellar image is finally achieved, with an order of mag- 
nitude reduction in any signature of fibre-to-fibre cross 
contamination. 

10 Conclusion 

We have presented an optimal extraction methodol- 
ogy for use with multi-fibre spectroscopy in the regime 
where fibres are tightly packed onto the detector array. 
We demonstrate that high accuracy can be achieved 
with this approach, minimising the impact of fibre-to- 
fibre cross-contamination. Additionally we show that 
the Gaussian least squares fitting extraction used by 
default for AAOmega-MOS data reduce within 2dfdr 
is adequate for the well-separated fibre profiles of this 
instrument mode provided the range of input target 
magnitudes is kept below Amag<3. 
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Figure 7: A reconstructed 2D image, created by collapsing the spectral data cube over wavelength, is 
shown for a SPIRAL-IFU observation of a standard star. Four dead fibres have been interpolated in the 
data, and are marked with small spots. 

Left - The reconstructed image obtained using a simple tramline extraction model. The vertical structure 
is not a telescope diffraction spike, but rather an artifact of fibre-to-fibre crosstalk. 

Centre - The Gaussian summation extraction suppresses noise but not the fibre-to-fibre cross contamina- 
tion. 

Right - The reconstructed image obtained using the optimal extraction model. The vertical extraction 
artifact structure removed and a rounded PSF is recovered. 
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